#ifndef __POLYNOMIAL_H__
#define __POLYNOMIAL_H__

#include "Math.h"

template <class Real>
class Polynomial
{
public:
	// construction and destruction
	Polynomial( int degree=-1 );
	Polynomial( const Polynomial& poly );
	~Polynomial();


	// member access
	void set_degree( int degree );
	int get_degree() const;
	operator const Real* () const;
	operator Real* ();

	Real operator[] (int i) const;
	Real& operator[] ( int i );

	// assignment
	Polynomial& operator= ( const Polynomial& poly );

	// evaluation
	Real operator() ( Real t ) const;

	// arithemtic operations
	Polynomial operator+ ( const Polynomial& poly ) const;
	Polynomial operator- ( const Polynomial& poly ) const;
	Polynomial operator* ( const Polynomial& poly ) const;
	Polynomial operator+ ( Real scalar ) const;
	Polynomial operator- ( Real scalar ) const;
	Polynomial operator* ( Real scalar ) const;
	Polynomial operator/ ( Real scalar ) const;
	Polynomial operator- ( );

	// arithmetic updates
	Polynomial& operator+= ( const Polynomial& poly );
	Polynomial& operator-= ( const Polynomial& poly );
	Polynomial& operator*= ( const Polynomial& poly );
	Polynomial& operator+= ( Real scalar );
	Polynomial& operator-= ( Real scalar );
	Polynomial& operator*= ( Real scalar );
	Polynomial& operator/= ( Real scalar );
	
	// derivative
	Polynomial derivative ( ) const;

	// inversion ( invpoly[i] = poly[degree-i] for 0 <= i <= degree )
	Polynomial inversion ( ) const;

	// Reduce degree by eliminating all (nearly) zero leading coefficients 
	// and by making the leading coefficient one. The input parameter is 
	// the threshold for specifiying that a coefficient is effectively zero.
	void compress( Real epsilon );


	// Divide P(t) with D(t) if degree(P) >= degree(D), then P(t) = Q(t)*D(t) + R(t)
	// where Q(t) is quotient with degree(Q) = degree(P) - degree(D) and R(t) is
	// remainder with degree(R) < degree(D). If routine is called with degree(P) < degree(D),
	// Q = 0 and R = P are returned. The value of epsilon is used as a threshold on coefficients
	// of the remainder polynomial. If smaller coefficient is assumed to be zero.
	void divide( const Polynomial& div, Polynomial& quot, Polynomial& rem, Real epsilon ) const;

protected:
	int degree;
	Real* coeff;
};

template <class Real>
Polynomial<Real> operator* ( Real scalar, const Polynomial<Real>& poly );
typedef Polynomial<float> Polynomialf;
typedef Polynomial<double> Polynomiald;

#endif